clear;
clc;
%warning off;
[Ypos, Yneg, Yzero,branchi,branchj,branches,Xp,X0,nodes,YNd_index]=Read(CSEE22);
location=3;
type=1;
%  输入短路节点号和短路类型，计算故障点相关参数
% location=input('请输入短路节点号：');
% type=input('请输入短路类型： 1：单相接地短路 2：双相短路 3：三相短路 4：两相接地短路 ');

Zp=GetZf(nodes,location,Ypos);
Zn=GetZf(nodes,location,Yneg);
Z0=GetZf(nodes,location,Yzero);

if type==1 %单相接地短路
    Ip=1/(Zp(location)+Zn(location)+Z0(location));
    In=Ip;
    If0=Ip;
elseif type==2 %双相短路
    Ip=1/(Zp(location)+Zn(location));
    In=Ip*(-1);
    If0=0;
elseif type==3 %三相短路
    Ip=1/Zp(location);
    In=0;
    If0=0;
elseif type==4 %两相接地短路
    tmp=Zn(location)+Z0(location);
    Ip=1/(Zp(location)+Zn(location)*Z0(location)/tmp);
    In=-Ip*Z0(location)/tmp;
    If0=-Ip*Zn(location)/tmp;
end
U=ones(nodes,1);

% 计算短路实际节点电压
Ui1=U-Zp*Ip;
Ui2=-Zn*In;
Ui0=-Z0*If0;

% 计算短路支路电流
Iij1=(Ui1(branchi)-Ui1(branchj))./Xp;
Iij2=(Ui2(branchi)-Ui2(branchj))./Xp;
Iij0=(Ui0(branchi)-Ui0(branchj))./X0;

%移相
% ------------------------------------------------------------------------------
% 考虑YNd变压器的相位偏移
% 由于不知道连接组别号，按YNd11连接组处理
k=ones(branches,1);

% 数据中低压侧（三角接）在后
Ui1(branchj(YNd_index))=Ui1(branchj(YNd_index)).*exp(1i*1/6*pi*k(YNd_index));
Ui2(branchj(YNd_index))=Ui2(branchj(YNd_index)).*exp(-1i*1/6*pi*k(YNd_index));


Iij1c=-1*Iij1(YNd_index).*exp(1i*1/6*pi*k(YNd_index));
Iij1=[Iij1;Iij1c];

Iij2c=-1*Iij2(YNd_index).*exp(-1i*1/6*pi*k(YNd_index));
Iij2=[Iij2;Iij2c];

Iij0=[Iij0;Iij0(YNd_index)];

%追加支路
branchi=[branchi;branchj(YNd_index)];
branchj=[branchj;branchi(YNd_index)];

%计算短路支路电流
a=exp(1i*2/3*pi);
Iija = Iij1 + Iij2+Iij0;
Iijb = Iij1*a^2 + Iij2*a+Iij0;
Iijc = Iij1*a + Iij2*a^2+Iij0;

% 计算短路节点电压
Ua=Ui1+Ui2+Ui0;
Ub=Ui1*a^2+Ui2*a+Ui0;
Uc=Ui1*a+Ui2*a^2+Ui0;


% 显示短路节点电流和电压
showout(Ip,In,If0,type,location,Ui1,Ui2,Ui0,nodes,branchi,branchj,Ua,Ub,Uc,Iija,Iijb,Iijc,Iij1,Iij2,Iij0);


%相关函数
%---------------------------------------------------------------------------------------------------------
%---------------------------------------------------------------------------------------------------------
function [Ypos, Yneg, Yzero,branchi,branchj,branches,Xp,X0,nodes,YNd_index] = Read(filename)
% 加载数据
%读入一个.mat文件，文件名为filename，文件包含两个矩阵，一个是线路数据，一个是发电机数据
%线路矩阵满足支路号	 支路类型(S) 	节点(I)	 节点(J) 	正序电抗(X)	零序电抗(Xo)
%式中支路类型：1：线路，2：YN,yn变压器，3：YN，d变压器
%发电机矩阵满足发电机节点(I) 	直轴次暂态电抗(xd") 	负序电抗(x2)
%全都为标幺值
%输出为三序节点导纳矩阵，支路起始节点和终止节点列表，支路数，节点数，正序阻抗和零序阻抗


linedata = filename.line;
gendata=filename.gen;

% 获取节点数和支路数
branch =  size(linedata);
branches = branch(1);
nodes = max(linedata(:,[2 3]),[],"all");

%提取出正序阻抗和零序阻抗
Xp=linedata(:,5)*1j;
X0=linedata(:,6)*1j;

% 提取出支路的起始节点和终止节点列表
branchi=linedata(:,3);
branchj=linedata(:,4);

% 计算线路的正负序节点导纳矩阵
positiveimpedance=linedata(:,5)*1j;
Ys = sparse(branchi, branchj, -1./positiveimpedance, nodes, nodes);
Ys = Ys+Ys.';

% 追加对角元素
Ys=Ys+sparse(branchi,branchi,1./positiveimpedance,nodes,nodes)+sparse(branchj,branchj,1./positiveimpedance,nodes,nodes);


% 计算添加发电机直轴电抗的正序节点导纳矩阵
rowi=gendata(:,1);
genx=sparse(rowi,rowi,1./(gendata(:,2)*1j),nodes,nodes);
Ypos=Ys+genx;

% 计算添加发电机负序电抗的负序节点导纳矩阵
gennx=sparse(rowi,rowi,1./(gendata(:,3)*1j),nodes,nodes);
Yneg=Ys+gennx;


% 计算零序节点导纳矩阵
YCom_index=[find(linedata(:,2)==1);find(linedata(:,2)==2)];%获得非Ynd变压器的支路索引
YNd_index=find(linedata(:,2)==3);%获得Ynd变压器的支路索引
zeroimpedance=linedata(YCom_index,6)*1j;
Yzero = sparse(branchi(YCom_index), branchj(YCom_index), -1./zeroimpedance, nodes, nodes);
Yzero = Yzero+Yzero.';  

Yzero=Yzero+sparse(branchi(YCom_index),branchi(YCom_index),1./zeroimpedance,nodes,nodes)+...
    sparse(branchj(YCom_index),branchj(YCom_index),1./zeroimpedance,nodes,nodes);
rowi=linedata(YNd_index,3);

  Yzero=Yzero+sparse(rowi,rowi,1./(linedata(YNd_index,6)*1j),nodes,nodes);
end

function Z=GetZf(nodes,location,Y)
% 计算故障阻抗阵的f列矩阵
StandI=sparse(location,1,1,nodes,1);
Z=Y\StandI;
Z(isnan(Z))=0;
end

function showout(Ip,In,If0,type,location,Ui1,Ui2,Ui0,nodes,branchi,branchj,Ua,Ub,Uc,Iija,Iijb,Iijc,Iij1,Iij2,Iij0)
%输出最终结果
% 显示一个标头，然后显示短路节点电流和电压
%读入短路点三序电流，短路类型，位置，三序电压，节点数，支路起始节点和终止节点列表，三相电压，三相电流
rowt={'单相接地短路' ,'两相短路','三相短路','两相接地'};
fprintf('----------%d节点%s电流和电压-----------\n',  location,string(rowt{type}));

colt={'故障类型  ' ,'  正序 '	,'  负序 '	,'  零序  ','  正序  ',	'  负序  '	,'   零序'};
format("short");
for index = 1:length(colt)
    fprintf('%-18s\t',colt{index});
end
fprintf('\n');
fprintf("%-10s",rowt{type});
out1=[Ip,In,If0,Ui1(location),Ui2(location),Ui0(location)];
out1=num2str(out1(1,:),'%15.4f');
disp(out1);
fprintf('\n-------------------------------------------------------------------------------------------------\n');

fprintf('\n----------短路各节点电压-----------\n');
fprintf('节点i\t  正序\t            负序\t           零序\t            A相\t               B相\t                C相\n');
U=full([Ui1,Ui2,Ui0 Ua,Ub,Uc]);

for index = 1:nodes
    fprintf('%-10d',index);
    Umodele=abs(U(index,:));
    Uangle=rad2deg(angle(U(index,:))) ;
    c=sprintf("%.3f∠%.2f°\t\t\t",[Umodele;Uangle]);
    fprintf(c+'\n');
    %disp(num2str(U(index,:),'%15.3f'));
end
fprintf('\n');
fprintf('\n----------短路各支路电流-----------\n');
fprintf('%-6s\t%-6s\t%-14s\t%-14s\t%-14s\t%-14s\t%-14s\t%-18s\n','支路i','支路j','正序电流','负序电流','零序电流','a相电流','b相电流','c相电流');

Iout=full([Iij1,Iij2,Iij0,Iija,Iijb,Iijc]);
for index = 1:length(branchi)
    fprintf('%-10d%-10d',branchi(index),branchj(index));
%     Imodele=abs(Iout(index,:));
%     Iangle=rad2deg(angle(Iout(index,:))) ;
%     c=sprintf("%.3f∠%.2f°\t\t\t",[Imodele;Iangle]);
    %fprintf(c+'\n');
    disp(num2str(Iout(index,:),'%15.3f'));
end
end